#include "fermiqcd.h"

int main(int argc, char** argv) {
  define_base_matrices("FERMILAB");
  mdp.open_wormholes(argc,argv);    
  int nc=3;
  int L[]={12,12,12,12};
  mdp_lattice lattice(4,L,default_partitioning0,torus_topology,0,2,false);
  gauge_field U(lattice,nc);
  char filename[128];
  mdp_site x(lattice);
  vector<mdp_real> p(4);
  for(int mu=0; mu<4; mu++) p[mu]=L[mu]/2;
  mdp_field<mdp_real> pion(lattice);
  pion = 0;

  // set parameters to naive quark propagator
  coefficients quark;
  quark["kappa_s"]=0.123;
  quark["kappa_t"]=0.123;
  quark["r_s"]=0.0;
  quark["r_t"]=0.0;
  quark["c_{sw}"]=0.0;
  default_fermi_action=FermiCloverActionFast::mul_Q;
  mdp_real absolute_precision = 10E-12;
  mdp_real relative_precision = 10E-8;
  
  // build a gauge config with a single instanton
  Instanton4D A(3,0,1, 1.0,3,p);      
  // set_cold(U);
  forallsites(x)
    for(int mu=0; mu<U.ndim; mu++)
      U(x,mu)=exp(-I*A(x,mu));
  
  // build source
  fermi_field chi(lattice,nc);
  fermi_field psi(lattice,nc);
  for(int beta=0; beta<4; beta++)
    for(int j=0; j<nc; j++) {
      forallsites(x)
	if(x(0)==p[0] && x(1)==p[1] && x(2)==p[2] && x(3)==p[3])
	  for(int alpha=0; alpha<4; alpha++)
	    for(int i=0; i<nc; i++)
	      psi(x,alpha,i)=(alpha==beta && i==j)?1:0;
	else
	  for(int alpha=0; alpha<4; alpha++)
	    for(int i=0; i<nc; i++)
	      psi(x,alpha,i)=0;
      psi.update();

      // compute em field if clover
      if(quark["c_{sw}"]!=0) compute_em_field(U);

      // invert Dirac action
      mul_invQ(chi,psi,U,quark,absolute_precision,relative_precision);        	

      // compute pion
      forallsites(x)
	for(int alpha=0; alpha<4; alpha++)
	  for(int i=0; i<nc; i++)
	    pion(x)+=pow(abs(chi(x,alpha,i)),2);
    }
  
  pion.save_vtk("pion_field.vtk");
  
  mdp.close_wormholes();
  return 0;
}
